tic
clear;

load 'C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM_R&R\Monte_Code\Results_Final\sp_tscs_phi50_rho25'; 

clear std; 

time = 10;

unit_1 = 51; 
unit_2 = 72; 

X = store_x0(nt-(n-1):nt,1);
X_prod = ones(t,1); 
X_steady = X*X_prod';
X_steady = reshape(X_steady,nt,1);

trials = 1; 

for h=1:trials

mult_fx = inv(i_nt - 0.25*WNT - 0.5*TL);
SIGMA = mult_fx*mult_fx';
rf_cut_full = mult_fx*(-1.5 + 3*X_steady); 

p_1 = zeros(time,1);
p0_1 = normcdf(rf_cut_full(unit_2),0,sqrt(SIGMA(unit_2,unit_2)));

for i=1:time
        p_1(i) = mvncdf([-Inf,-Inf], [rf_cut_full([unit_2]), rf_cut_full([n*(i-1) + unit_1])],0,SIGMA([unit_2,n*(i-1) + unit_1],[unit_2,n*(i-1) + unit_1])); 
end

cp_1 = zeros(time,1);


for i = 1:time
    cp_1(i) = p_1(i)/p0_1;
end

p_0 = zeros(time,1);
p0_0 = normcdf(-rf_cut_full(unit_2),0,sqrt(SIGMA(unit_2,unit_2)));

for i=1:time
        p_0(i) = mvncdf([rf_cut_full([unit_2]),-Inf], [Inf, rf_cut_full([n*(i-1) + unit_1])],0,SIGMA([unit_2,n*(i-1) + unit_1],[unit_2,n*(i-1) + unit_1])); 
end

cp_0 = zeros(time,1);


for i = 1:time
    cp_0(i) = p_0(i)/p0_0;
end
cp_diff = cp_1-cp_0;

cp_diff_true(:,h) = cp_diff;

end

true_fx = cp_diff_true;

save 'C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM_R&R\Monte_Code\Results_Final\sp_tscs_phi50_rho25_fx_true'; 